** Replication files for Weisiger and Gartzke "Debating the Democratic Peace in 
** the International System"

** We work from the replication files provided by Crescenzi and Kadera, which in 
** turn are based on our own replication files from GW 2014.  The final figure
** is constructed using the replication data provided by Choi.

** This do file replicates both core models and extensions, as follows:
** 1) Replication of our tables 1, 2, and 3
** 2) Substantiating the claim that Democratic Community is negative and 
**    statistically significant in half of models omitting Systemic Difference
**    but in no model that includes Difference.
** 3) Replicating figure 1
** 4) Replicating figure 2

use crescenzi-kadera-2015-ISQ-Table2.dta, clear
set more off

** Replicating Table 1

** Their table 1, model C
nbreg fat1 dempower dempower2 diff1 pcenerg sysdepstate states year, robust
eststo ck1

** Dropping trade to recover the full set of observations
nbreg fat1 dempower dempower2 diff1 pcenerg states year, robust
eststo ck2

** Dropping the quadratic term
nbreg fat1 dempower diff1 pcenerg sysdepstate states year, robust
eststo ck3

lab var dempower "Democratic Power"
lab var dempower2 "Democratic Power$^2$"
lab var diff1 "Difference"
lab var pcenerg "Development"
lab var sysdepstate "Trade"
lab var states "\# of States"
lab var year "Year"

esttab ck1 ck2 ck3 using ckresults.tex, tex b(a2) se l mtitles("Model 1" "Model 2" "Model 3") title("Replications of Crescenzi and Kadera (Table 2)") nonotes addnotes("\sym{*} \(p<0.05\), \sym{**} \(p<0.01\)") replace

** Replicating Table 2

nbreg fat1 regstrength pcenerg states year, nolog robust
eststo ck41
nbreg fat1 polave pcenerg states year, nolog robust
eststo ck42
nbreg fat1 regstrength diff1 pcenerg states year, nolog robust
eststo ck43

nbreg fat1 c_regstrength c_regstrsq pcenerg states year, nolog robust
eststo ck51
nbreg fat1 polave polave2 pcenerg states year, nolog robust
eststo ck52
nbreg fat1 c_regstrength c_regstrsq diff1 pcenerg states year, nolog robust
eststo ck53

nbreg fat1 regstrength pcenerg sysdepstate states year, nolog robust
eststo ck61
nbreg fat1 polave pcenerg sysdepstate states year, nolog robust
eststo ck62
nbreg fat1 regstrength diff1 pcenerg sysdepstate states year, nolog robust
eststo ck63

lab var regstrength "Democratic Community"
lab var polave "Average Polity"
lab var polave2 "Average Polity$^2$"
lab var diff1 "Difference"
lab var pcenerg "Development"
lab var sysdepstate "Trade"
lab var c_regstrsq "Democratic Community$^2$"
lab var states "\# of States"
lab var year "Year"

esttab ck41 ck42 ck43 ck51 ck52 ck53 ck61 ck62 ck63 using ckresults.tex, tex b(a2) se l mtitles("Model 1" "Model 2" "Model 3") title("Replications of Crescenzi and Kadera (Table 2)") nonotes addnotes("\sym{*} \(p<0.05\), \sym{**} \(p<0.01\)") append

** Replicating Table 3

use crescenzi-kadera-2015-ISQ-Table3, clear

** First replicating Crescenzi and Kadera table 3, model 3

logit deadlyl c_regstrength c_regstrengthsq pcenerg demloi engypop logdist cntgdumy allydumy capratio onemajor deadyrs deadyer*, cluster(dyadid)
eststo tab31

** Now restoring systemic difference

logit deadlyl c_regstrength c_regstrengthsq diff1 pcenerg demloi engypop logdist cntgdumy allydumy capratio onemajor deadyrs deadyer*, cluster(dyadid)
eststo tab32

lab var c_regstrength "Democratic Community"
lab var c_regstrengthsq "Democratic Community$^2$"
lab var diff1 "Systemic Difference"
lab var pcenerg "Systemic Development"
lab var demloi "Democracy (low)"
lab var engypop "Dyadic Development"
lab var logdist "Distance"
lab var cntgdumy "Contiguity"
lab var allydumy "Alliance"
lab var capratio "Relative Capabilities"
lab var onemajor "Major Power"

esttab tab31 tab32 using ckresults.tex, tex b(a2) se l mtitles("Model 1" "Model 2") title("Replication of Crescenzi and Kadera (Table 3)") nonotes addnotes("\sym{*} \(p<0.05\), \sym{**} \(p<0.01\)") d(dead*) o(c_* diff1) append

** The remainder of the do file contains code to generate the figures included
** in the response and substantiates a couple specific claims that we make.

** We note in the response that Democratic Community has a negative and 
** significant relationship with systemic conflict levels in 8 of 16 
** specifications when omitting Systemic Difference, but is never significant
** when Difference is included in the regression.  The next section validates
** that claim.  The loop is designed to (marginally) abbreviate the code 
** necessary to run the 32 regressions and store the relevant information.

use crescenzi-kadera-2015-ISQ-Table2.dta, clear

keep fat1 pcenerg wgdppc3 states year diff1 sysdepstate polave* propdem* dempower* regstrength

local devlist="pcenerg"
local controllist="diff1 states sysdepstate year"
local demlist="polave regstrength"

gen dev=""
gen dem=""
gen controls=""
foreach var in dev dem {
	gen `var'beta=.
	gen `var'hi=.
	gen `var'lo=.
}

local i=1
foreach dev in pcenerg {
*	di "Now moving to development variable `dev'"
	foreach dem in polave regstrength {
		di "Now moving to democracy variable `dem'"
		qui nbreg fat1 `dev' `dem', nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dev="`dev'" in `i'
		qui replace dem="`dem'" in `i'
		qui replace devbeta=A[1,1] in `i'
		qui replace devhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace devlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace dembeta=A[1,2] in `i'
		qui replace demhi=A[1,2]+invnormal(.975)*sqrt(B[2,2]) in `i'
		qui replace demlo=A[1,2]-invnormal(.975)*sqrt(B[2,2]) in `i'
		local i=`i'+1
		foreach control in diff1 states sysdepstate year {
			qui nbreg fat1 `dev' `dem' `control', nolog robust
			mat def A=e(b)
			mat def B=e(V)
			qui replace dev="`dev'" in `i'
			qui replace dem="`dem'" in `i'
			qui replace controls="`control'" in `i'
			qui replace devbeta=A[1,1] in `i'
			qui replace devhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
			qui replace devlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
			qui replace dembeta=A[1,2] in `i'
			qui replace demhi=A[1,2]+invnormal(.975)*sqrt(B[2,2]) in `i'
			qui replace demlo=A[1,2]-invnormal(.975)*sqrt(B[2,2]) in `i'
			local i=`i'+1	
		}
		foreach control in states sysdepstate year {
			qui nbreg fat1 `dev' `dem' diff1 `control', nolog robust
			mat def A=e(b)
			mat def B=e(V)
			qui replace dev="`dev'" in `i'
			qui replace dem="`dem'" in `i'
			qui replace controls="diff1 `control'" in `i'
			qui replace devbeta=A[1,1] in `i'
			qui replace devhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
			qui replace devlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
			qui replace dembeta=A[1,2] in `i'
			qui replace demhi=A[1,2]+invnormal(.975)*sqrt(B[2,2]) in `i'
			qui replace demlo=A[1,2]-invnormal(.975)*sqrt(B[2,2]) in `i'
			local i=`i'+1
		}
		foreach control in sysdepstate year {
			qui nbreg fat1 `dev' `dem' states `control', nolog robust
			mat def A=e(b)
			mat def B=e(V)
			qui replace dev="`dev'" in `i'
			qui replace dem="`dem'" in `i'
			qui replace controls="states `control'" in `i'
			qui replace devbeta=A[1,1] in `i'
			qui replace devhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
			qui replace devlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
			qui replace dembeta=A[1,2] in `i'
			qui replace demhi=A[1,2]+invnormal(.975)*sqrt(B[2,2]) in `i'
			qui replace demlo=A[1,2]-invnormal(.975)*sqrt(B[2,2]) in `i'
			local i=`i'+1
		}
		qui nbreg fat1 `dev' `dem' sysdepstate year, nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dev="`dev'" in `i'
		qui replace dem="`dem'" in `i'
		qui replace controls="sysdepstate year" in `i'
		qui replace devbeta=A[1,1] in `i'
		qui replace devhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace devlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace dembeta=A[1,2] in `i'
		qui replace demhi=A[1,2]+invnormal(.975)*sqrt(B[2,2]) in `i'
		qui replace demlo=A[1,2]-invnormal(.975)*sqrt(B[2,2]) in `i'
		local i=`i'+1
		qui nbreg fat1 `dev' `dem' diff1 states sysdepstate, nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dev="`dev'" in `i'
		qui replace dem="`dem'" in `i'
		qui replace controls="diff1 states sysdepstate" in `i'
		qui replace devbeta=A[1,1] in `i'
		qui replace devhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace devlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace dembeta=A[1,2] in `i'
		qui replace demhi=A[1,2]+invnormal(.975)*sqrt(B[2,2]) in `i'
		qui replace demlo=A[1,2]-invnormal(.975)*sqrt(B[2,2]) in `i'
		local i=`i'+1
		qui nbreg fat1 `dev' `dem' diff1 states year, nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dev="`dev'" in `i'
		qui replace dem="`dem'" in `i'
		qui replace controls="diff1 states year" in `i'
		qui replace devbeta=A[1,1] in `i'
		qui replace devhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace devlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace dembeta=A[1,2] in `i'
		qui replace demhi=A[1,2]+invnormal(.975)*sqrt(B[2,2]) in `i'
		qui replace demlo=A[1,2]-invnormal(.975)*sqrt(B[2,2]) in `i'
		local i=`i'+1
		if "`dem'"=="dempower" & "`dev'"=="wgdppc3" { /* does not converge */
		}
		else {
		qui nbreg fat1 `dev' `dem' diff1 sysdepstate year, nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dev="`dev'" in `i'
		qui replace dem="`dem'" in `i'
		qui replace controls="diff1 sysdepstate year" in `i'
		qui replace devbeta=A[1,1] in `i'
		qui replace devhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace devlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace dembeta=A[1,2] in `i'
		qui replace demhi=A[1,2]+invnormal(.975)*sqrt(B[2,2]) in `i'
		qui replace demlo=A[1,2]-invnormal(.975)*sqrt(B[2,2]) in `i'
		local i=`i'+1
		}
		qui nbreg fat1 `dev' `dem' states sysdepstate year, nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dev="`dev'" in `i'
		qui replace dem="`dem'" in `i'
		qui replace controls="states sysdepstate year" in `i'
		qui replace devbeta=A[1,1] in `i'
		qui replace devhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace devlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace dembeta=A[1,2] in `i'
		qui replace demhi=A[1,2]+invnormal(.975)*sqrt(B[2,2]) in `i'
		qui replace demlo=A[1,2]-invnormal(.975)*sqrt(B[2,2]) in `i'
		local i=`i'+1
		qui nbreg fat1 `dev' `dem' diff1 states sysdepstate year, nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dev="`dev'" in `i'
		qui replace dem="`dem'" in `i'
		qui replace controls="diff1 states sysdepstate year" in `i'
		qui replace devbeta=A[1,1] in `i'
		qui replace devhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace devlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace dembeta=A[1,2] in `i'
		qui replace demhi=A[1,2]+invnormal(.975)*sqrt(B[2,2]) in `i'
		qui replace demlo=A[1,2]-invnormal(.975)*sqrt(B[2,2]) in `i'
		local i=`i'+1
	}
	** And now the same set of regressions, but without pcenerg
		foreach dem in polave regstrength {
		di "Now moving to democracy variable `dem'"
		qui nbreg fat1 `dem', nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dem="`dem'" in `i'
		qui replace dembeta=A[1,1] in `i'
		qui replace demhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace demlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		local i=`i'+1
		foreach control in diff1 states sysdepstate year {
			qui nbreg fat1 `dem' `control', nolog robust
			mat def A=e(b)
			mat def B=e(V)
			qui replace dem="`dem'" in `i'
			qui replace controls="`control'" in `i'
			qui replace dembeta=A[1,1] in `i'
			qui replace demhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
			qui replace demlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
			local i=`i'+1	
		}
		foreach control in states sysdepstate year {
			qui nbreg fat1 `dem' diff1 `control', nolog robust
			mat def A=e(b)
			mat def B=e(V)
			qui replace dem="`dem'" in `i'
			qui replace controls="diff1 `control'" in `i'
			qui replace dembeta=A[1,1] in `i'
			qui replace demhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
			qui replace demlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
			local i=`i'+1
		}
		foreach control in sysdepstate year {
			qui nbreg fat1 `dev' `dem' states `control', nolog robust
			mat def A=e(b)
			mat def B=e(V)
			qui replace dem="`dem'" in `i'
			qui replace controls="states `control'" in `i'
			qui replace dembeta=A[1,1] in `i'
			qui replace demhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
			qui replace demlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
			local i=`i'+1
		}
		qui nbreg fat1 `dem' sysdepstate year, nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dem="`dem'" in `i'
		qui replace controls="sysdepstate year" in `i'
		qui replace dembeta=A[1,1] in `i'
		qui replace demhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace demlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		local i=`i'+1
		qui nbreg fat1 `dem' diff1 states sysdepstate, nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dem="`dem'" in `i'
		qui replace controls="diff1 states sysdepstate" in `i'
		qui replace dembeta=A[1,1] in `i'
		qui replace demhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace demlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		local i=`i'+1
		qui nbreg fat1 `dem' diff1 states year, nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dem="`dem'" in `i'
		qui replace controls="diff1 states year" in `i'
		qui replace dembeta=A[1,1] in `i'
		qui replace demhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace demlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		local i=`i'+1
		qui nbreg fat1 `dem' diff1 sysdepstate year, nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dem="`dem'" in `i'
		qui replace controls="diff1 sysdepstate year" in `i'
		qui replace dembeta=A[1,1] in `i'
		qui replace demhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace demlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		local i=`i'+1
		qui nbreg fat1 `dem' states sysdepstate year, nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dem="`dem'" in `i'
		qui replace controls="states sysdepstate year" in `i'
		qui replace dembeta=A[1,1] in `i'
		qui replace demhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace demlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		local i=`i'+1
		qui nbreg fat1 `dem' diff1 states sysdepstate year, nolog robust
		mat def A=e(b)
		mat def B=e(V)
		qui replace dem="`dem'" in `i'
		qui replace controls="diff1 states sysdepstate year" in `i'
		qui replace dembeta=A[1,1] in `i'
		qui replace demhi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
		qui replace demlo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
		local i=`i'+1
	}
}

** Total number of regressions involving regstrength (Democratic Community)
count if dem=="regstrength"
** Total number of regressions involving regstrength but not systemic difference
count if dem=="regstrength" & strpos(controls,"diff1")==0
** Regressions in which regstrength is negative and significant, when systemic
** difference is not included as a control
count if dem=="regstrength" & demhi<0 & strpos(controls,"diff1")==0
** Regressions in which regstrength is negative and significant, when systemic
** difference is included as a control
count if dem=="regstrength" & demhi<0 & strpos(controls,"diff1")>0


** Generating the data used to create figure 1.  This section is commented out 
** because it takes a long time to run -- if run, it will produce ck_preds.dta.

/*

use crescenzi-kadera-2015-ISQ-Table3.dta, clear
set more off

gen n=_n
gen yr=1815+_n if _n<=185
foreach var in polave regstrength {
	foreach dyad in bothdem bothaut mixed {
		foreach form in lin quad {
			foreach ext in b hi lo {
				gen pred_`var'_`form'_`dyad'_`ext'=.
			}
		}
	}
}

logit deadlyl c_polave c_polavesq pcenerg diff1 demloi engypop dydiff logdist cntgdumy allydumy capratio onemajor deadyrs deadyer*, cluster(dyadid)
local i=1
forvalues yr=1816/2000 {
	if `yr'/25==round(`yr'/25) di `yr'
	qui sum n if year==`yr'
	local n=r(min)
	local polave=c_polave[`n']
	local polave2=c_polavesq[`n']
	local diff1=diff1[`n']
	qui margins, at((median) _all c_polave=`polave' c_polavesq=`polave2' diff1=`diff1' demloi=10 dydiff=0)
	mat def A=r(b)
	mat def B=r(V)
	qui replace pred_polave_quad_bothdem_b=A[1,1] in `i'
	qui replace pred_polave_quad_bothdem_hi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
	qui replace pred_polave_quad_bothdem_lo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
	qui margins, at((median) _all c_polave=`polave' c_polavesq=`polave2' diff1=`diff1' demloi=0 dydiff=0)
	mat def A=r(b)
	mat def B=r(V)
	qui replace pred_polave_quad_bothaut_b=A[1,1] in `i'
	qui replace pred_polave_quad_bothaut_hi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
	qui replace pred_polave_quad_bothaut_lo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
	qui margins, at((median) _all c_polave=`polave' c_polavesq=`polave2' diff1=`diff1' demloi=0 dydiff=10)
	mat def A=r(b)
	mat def B=r(V)
	qui replace pred_polave_quad_mixed_b=A[1,1] in `i'
	qui replace pred_polave_quad_mixed_hi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
	qui replace pred_polave_quad_mixed_lo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
	local i=`i'+1
}
save ck_preds, replace

logit deadlyl c_polave pcenerg diff1 demloi engypop dydiff logdist cntgdumy allydumy capratio onemajor deadyrs deadyer*, cluster(dyadid)
local i=1
forvalues yr=1816/2000 {
	if `yr'/25==round(`yr'/25) di `yr'
	qui sum n if year==`yr'
	local n=r(min)
	local polave=c_polave[`n']
	local diff1=diff1[`n']
	qui margins, at((median) _all c_polave=`polave' diff1=`diff1' demloi=10 dydiff=0)
	mat def A=r(b)
	mat def B=r(V)
	qui replace pred_polave_lin_bothdem_b=A[1,1] in `i'
	qui replace pred_polave_lin_bothdem_hi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
	qui replace pred_polave_lin_bothdem_lo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
	qui margins, at((median) _all c_polave=`polave' diff1=`diff1' demloi=0 dydiff=0)
	mat def A=r(b)
	mat def B=r(V)
	qui replace pred_polave_lin_bothaut_b=A[1,1] in `i'
	qui replace pred_polave_lin_bothaut_hi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
	qui replace pred_polave_lin_bothaut_lo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
	qui margins, at((median) _all c_polave=`polave' diff1=`diff1' demloi=0 dydiff=10)
	mat def A=r(b)
	mat def B=r(V)
	qui replace pred_polave_lin_mixed_b=A[1,1] in `i'
	qui replace pred_polave_lin_mixed_hi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
	qui replace pred_polave_lin_mixed_lo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
	local i=`i'+1
}
save ck_preds, replace

logit deadlyl c_regstrength c_regstrengthsq pcenerg diff1 demloi engypop dydiff logdist cntgdumy allydumy capratio onemajor deadyrs deadyer*, cluster(dyadid)
local i=1
forvalues yr=1816/2000 {
	if `yr'/25==round(`yr'/25) di `yr'
	qui sum n if year==`yr'
	local n=r(min)
	local polave=c_polave[`n']
	local polave2=c_polavesq[`n']
	local diff1=diff1[`n']
	qui margins, at((median) _all c_regstrength=`polave' c_regstrengthsq=`polave2' diff1=`diff1' demloi=10 dydiff=0)
	mat def A=r(b)
	mat def B=r(V)
	qui replace pred_regstrength_quad_bothdem_b=A[1,1] in `i'
	qui replace pred_regstrength_quad_bothdem_hi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
	qui replace pred_regstrength_quad_bothdem_lo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
	qui margins, at((median) _all c_regstrength=`polave' c_regstrengthsq=`polave2' diff1=`diff1' demloi=0 dydiff=0)
	mat def A=r(b)
	mat def B=r(V)
	qui replace pred_regstrength_quad_bothaut_b=A[1,1] in `i'
	qui replace pred_regstrength_quad_bothaut_hi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
	qui replace pred_regstrength_quad_bothaut_lo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
	qui margins, at((median) _all c_regstrength=`polave' c_regstrengthsq=`polave2' diff1=`diff1' demloi=0 dydiff=10)
	mat def A=r(b)
	mat def B=r(V)
	qui replace pred_regstrength_quad_mixed_b=A[1,1] in `i'
	qui replace pred_regstrength_quad_mixed_hi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
	qui replace pred_regstrength_quad_mixed_lo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
	local i=`i'+1
}
save ck_preds, replace

logit deadlyl c_regstrength pcenerg diff1 demloi engypop dydiff logdist cntgdumy allydumy capratio onemajor deadyrs deadyer*, cluster(dyadid)
local i=1
forvalues yr=1816/2000 {
	if `yr'/25==round(`yr'/25) di `yr'
	qui sum n if year==`yr'
	local n=r(min)
	local polave=c_polave[`n']
	local diff1=diff1[`n']
	qui margins, at((median) _all c_regstrength=`polave' diff1=`diff1' demloi=10 dydiff=0)
	mat def A=r(b)
	mat def B=r(V)
	qui replace pred_regstrength_lin_bothdem_b=A[1,1] in `i'
	qui replace pred_regstrength_lin_bothdem_hi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
	qui replace pred_regstrength_lin_bothdem_lo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
	qui margins, at((median) _all c_regstrength=`polave' diff1=`diff1' demloi=0 dydiff=0)
	mat def A=r(b)
	mat def B=r(V)
	qui replace pred_regstrength_lin_bothaut_b=A[1,1] in `i'
	qui replace pred_regstrength_lin_bothaut_hi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
	qui replace pred_regstrength_lin_bothaut_lo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
	qui margins, at((median) _all c_regstrength=`polave' diff1=`diff1' demloi=0 dydiff=10)
	mat def A=r(b)
	mat def B=r(V)
	qui replace pred_regstrength_lin_mixed_b=A[1,1] in `i'
	qui replace pred_regstrength_lin_mixed_hi=A[1,1]+invnormal(.975)*sqrt(B[1,1]) in `i'
	qui replace pred_regstrength_lin_mixed_lo=A[1,1]-invnormal(.975)*sqrt(B[1,1]) in `i'
	local i=`i'+1
}

save ck_preds, replace

use ck_preds.dta, clear

#del ;
twoway rarea pred_polave_lin_mixed_hi pred_polave_lin_mixed_lo yr, color(gs14) ||
rarea pred_polave_lin_bothaut_hi pred_polave_lin_bothaut_lo yr, color(gs14) || 
line pred_polave_lin_mixed_b pred_polave_lin_bothaut_b yr, lc(black black)
graphregion(fcolor(white)) legend(off) title("Average Polity" "Linear Term Only")
xtitle("Year") text(.0046 1945 "Mixed Dyad", place(e)) 
text(.0003 1925 "Autocratic Dyad", place(e)) name(polave_lin, replace)
ylabel(.001(.002).005)
;

twoway rarea pred_polave_quad_mixed_hi pred_polave_quad_mixed_lo yr, color(gs14) ||
rarea pred_polave_quad_bothaut_hi pred_polave_quad_bothaut_lo yr, color(gs14) || 
line pred_polave_quad_mixed_b pred_polave_quad_bothaut_b yr, lc(black black)
graphregion(fcolor(white)) legend(off) title("Average Polity" "Linear and Quadratic")
xtitle("Year") text(.004 1945 "Mixed Dyad", place(e)) 
text(.0003 1925 "Autocratic Dyad", place(e)) name(polave_quad, replace)
ylabel(.001(.002).005)
;

twoway rarea pred_regstrength_quad_mixed_hi pred_regstrength_quad_mixed_lo yr, color(gs14) ||
rarea pred_regstrength_quad_bothaut_hi pred_regstrength_quad_bothaut_lo yr, color(gs14) || 
line pred_regstrength_quad_mixed_b pred_regstrength_quad_bothaut_b yr, lc(black black)
graphregion(fcolor(white)) legend(off) title("Democratic Community" "Linear and Quadratic")
xtitle("Year") text(.0045 1945 "Mixed Dyad", place(e)) 
text(.0003 1925 "Autocratic Dyad", place(e)) name(regstrength_quad, replace)
ylabel(.001(.002).005)
;

twoway rarea pred_regstrength_lin_mixed_hi pred_regstrength_lin_mixed_lo yr, color(gs14) ||
rarea pred_regstrength_lin_bothaut_hi pred_regstrength_lin_bothaut_lo yr, color(gs14) || 
line pred_regstrength_lin_mixed_b pred_regstrength_lin_bothaut_b yr, lc(black black)
graphregion(fcolor(white)) legend(off) title("Democratic Community" "Linear Term Only")
xtitle("Year") text(.0045 1945 "Mixed Dyad", place(e)) 
text(.0003 1925 "Autocratic Dyad", place(e)) name(regstrength_lin, replace)
ylabel(.001(.002).005)
;

graph combine regstrength_lin regstrength_quad polave_lin polave_quad, cols(2) 
graphregion(fcolor(white))
;

*/


** Generating figure 2.  Here we are using Choi's replication data.

use gartzke_weisiger_ISQ_2014_table_3_data_012014_skewness.dta, clear

des sysdems

#del ;
line sysdems diff1 year if dyadid==2200, lp(solid dash) graphregion(fcolor(white)) 
legend(label(1 "Skewness") label(2 "Standard Deviation")) xtitle("Year");

